Numerical modelling of elastic-plastic deformation at crack tips in composite material under stress wave loading

The paper deals with the elastic-plastic wave propagations in an anisotropic material subjected to impact loading with emphasis on an ideal fibre-reinforced cubic material under two-dimensional plane strain condition. The system of partial differential equations is hyperbolic. Therefore a characteristic-based difference method is proposed to calculate the formation of the plastic zone at the crack t ip. 1. I N T R O D U C T I O N It is well-known that numerous materials, e.g. woods and several metals behave anisotropic under elastic and plastic deformations, see e.g. [1] and [2]. The rapid development and use of fibrereinforced materials has essentially promoted the development of constitutive relations for anisotropic materials, see e.g. Rogers' review paper [3]. Those constitutive equations can be used in numerical computational analyses of stresses and strains. One of the important subjects in the analysis of stress is the dynamic behaviour of materials under strong impact loading conditions. The governing equations form a system of hyperbolic PDEs, and the effects of wave propagations cannot be neglected. It has been shown in [4,5] that the finite difference method based on the bicharacteristic analysis represents a good tool in solving those systems, since it can give a direct physical interpretation of the involved wave pattern. This paper extends the work presented in [5]. A finite difference scheme is formulated in order to model the stress wave propagation in a composite cubic material and to treat crack problems in such materials. Article published online by EDP Sciences and available at http://dx.doi.org/10.1051/jp4:1994807 03-54 JOURNAL DE PHYSIQUE IV 2. CONSTITUTIVE FORMULATION OF CUBIC MATERIALS As usual, the deformation of cubic materials can be divided into two parts: volume changes and shape changes. The volume changes can be assumed as elastic for most problems unless the solid is subjected extremely high pressure and temperature. Denoting by K the bulk modulus, the volume change may be expressed as 1 3K(ux + uy + a,) = EX + Ey + E,(1) The shape changes are usually divided into the elastic change and plastic change. In the incremental formulation of the theory of plasticity, F&) = d ( ~ : ~ B ~ S , ~ ) + dcfj, (2) where, E = (E, + + ~ , ) / 3 . The elastic deformation E; vanishes when the stresses are removed. For this case, Hooke's law holds, where p and G are the two shear moduli. In genera1 one has p # G for cubic material. The plastic deformation cTj is permanent and remains when the stress is removed. The increment d ~ ! ~ may be obtained using a flow rule with a scalar plastic potential F(uij), I 0, when F < 1 or F = l, d F < 0, dc$ = 1 i3F dx -when ac, 2 ' Suppose that the plastic potential F follows the von Mises' yield stress function obtained by Zheng and Betten [6] , where 1 12 = Z[(ux + (U' + (Ui o.)~], T2 = r,y + T& + r:=. The symbols Y and K represent the yield stresses for simple tension along a proper axis direction and simple shear transverse to the axis direction, respectively. In the isotropic case, Y = &c. But in the anisotropic case Y and K are independent. For composite materials with very strong fibres in the axis directions, Y may be much greater than K. In other words, only the value of T will cause plastic deformation, while the contribution from I can be neglected. With F from eq. ( 5 ) , the incremental components of plastic strain may be expressed by The multiplier dx relies on the hardening law of the material. Spencer [7] has proposed a concept of proportional hardening. According to this concept, the yield function F retains its form as a function of the stress invariants I and T during the plastic flow, but the parameters Y and n depend on the deformation history. In order to determine dx, the total differentiation is taken for eq. (5) along the yield surface F = 1, This suggests that the functional form for the yield stresses Y and K should be found first. An assumption has been made by Rogers in [3] for fibre-reinforced materials in which Y is affected only by E:, E; and efl, while n is affected only by 7:y, 7iz and yzz. If this assumption is applied to the present case, the functional forms for Y and K, will be For one-dimensional tension along the z-direction, i.e. a, = uy = 0, U, = Y, where E, = da,/ds, is the slope of the simple stress-strain curve in the plastic range. Applying this condition to eqs. ( 6 ) , and being aware of the change of plastic volume d(e; + &P, + efl) = 0, one gets de; in term of dx. Therefore, dY=E(Y)dX. Y (10) By the same way one gets dn considering a simple shear condition, where pp = d~,,/d&,, is the slope of simple shear curve in the plastic range. Therefore, d~ is obtained when dY and dn are substituted into eq. (7), Finally, the constitutive relations for the shape changes can be obtained combining eqs. (3) and (6); only three of those equations are written down here which will be used for the plane strain problem in the next section: l a -d(u, o,) + -(c, u , ) ~ x = d ( ~ E , ) , 2G 2Y2 C8-56 JOURNAL DE PHYSIQUE IV 3. THE PLANE STRAIN PROBLEM Suppose the deformation of a cubic material is restricted by E, = 0, . y, = 0 and .yy, = 0. The non-zero stress components are denoted by p = U,, q = U,, T = U, and T = T,,. The velocity components in zand y-directions are represented by u and v, and p is the mass density of the material. Then the governing equations become au ap a7 p a , = -+-, az ay A two-step numerical scheme solving the system (14) can be built up by the method introduced in [5]. In the first step, the flux in a grid point is calculated. Denote the cell (i, j ) by cell 1, cell (i, j + 1) by cell 2, cell (i + l, j ) by cell 3, cell (i + l, j + 1) by cell 4. The value of W in cell k at the time level tn is denoted by wk, and the unknown value in the grid point (i + $, j + $) at time level tn+) is simply denoted by W itself. Since the first three equations of eqs. (14) are linear, U , v and p + q + r can be calculated directly. But p T, q T and T need an iterative procedure. Denote P = p T, Q = q r, from last three equations of the system the formulas of the iteration becomes: where X = AtlAz, and the point (I,', T;), i.e. (p; T;, q; T;, 7;) lies on the yield function surface of the h-th cell, I,' 2 7; 2 F = (-) + (-1 =l . yk R k In the second step, the unknown functions in the cell's centers are updated. Again, the stress components p T, q T and T should be obtained by iteration. Figure 1: Plastic zone and yield stress contours Y/YO (solid lines) and n/no (dashed lines) for the semi-infinite crack problem of a cubic material with G/p = 1.2. The contour values are l(outer), 1.05, 1.1, 1.2 and 1.5(inner). 4. FORMATION OF PLASTIC ZONE AT A CRACK TIP The above described model has been used in order to calculate the plastic zone formation at a crack tip in a cubic material under plane strain. The crack is assumed as semi-infinite. The material parameters are chosen as p = 1, K = 7/15, p = 113 (which result in cl = 1 and c2 = l/&). The initial yield stress no = 1, and the plastic shear modulus (defined in eq. (11)) p(n) = 0.0625. The anisotropic properties are described by for the modelling of the discussed fibre-reinforced material. The computing domain is divided into 600 X 300 cells in Xand y-direction, respectively. The crack lies in the x-axis and takes 300 cells in length. The initial conditions are described by an incoming pair of tensional waves with which are symmetrical to the x-axis v q0 = v0 = 1.51554, u0 = = 0, = = l v !I0, (18) where, v = (3K 2G)/(6K + 2G) is Poisson's ratio. The amplitudes of the oncoming waves are Io = 0.35Y0. The contours of yield stress Y/YO and ~ / K O are plotted in Fig. 1. The results show that the elastic-plastic boundaries defined by the two yield stresses are approximately in agreement, but the contour lines inside the plastic zone are different. The orientation is also significant. The considered material is strengthened by fibres in X and y-directions, the plastic zone is stretched in one fibre direction. Being aware of sketched wave fronts in the fibre-reinforced material, one may recognize the fact that, although the longitudinal wave propagates in the fibre direction faster than in any other direction, the transverse wave (and therefore the slow plastic wave) becomes slower in the fibre's direction. For the inhomogeneous wave, the main part of energy is transported by the transverse wave and slow plastic wave. Therefore, the energy will be accumulated in the fibre direction which causes the plastic yield. This explanation can be confirmed considering another example, in which the two families of fibres are orientated under 45 and 135 degrees with z-axis, respectively. The material parameters JOURNAL D E PHYSIQUE 1V Figure 2: Plastic zone and yield stress contours Y/YO (solid lines) and &/no (dashed lines) for the semi-infinite crack problem of a cubic material with G/p = 0.3137. The contour values are l(outer), 1.05, 1.1, 1.2 and 1.5(inner). p, K, p, KO and p(%) are chosen as before, but the anisotropic parameters are chosen as The incoming wave amplitude is still controlled by I0 = 0.35Y0. In order to calculate this problem, the same program has been used with the CFL number clAt/Ax = 0.8. Although an overshut was caused in the elastic shock wave front, little influence could be recognized in the plastic zone where the unknown functions change smoothly. The results are shown in Fig. 2. This figure shows indeed that the plastic zone extends in the fibre directions.